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Abstract 

An algorithm is presented for unsteady two dimensional incompressible Navier-Stokes 
calculations. This algorithm is based on the fourth order Partial Differential Equation for 
incompressible fluid flow which uses the streamfunction as the only dependent variable. 
The algorithm is second order accurate in both time and space, it uses a multigrid solver at 
each time step, it is extremely efficient with respect to use of both CPU time and physical 
memory, and it is extremely robust with respect to Reynolds number. 
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1. SUMMARY 


An algorithm is presented for unsteady two dimensional incompressible Navier-Stokes 
calculations. This algorithm is based on the fourth order Partial Differential Equation for 
incompressible fluid flow which uses the streamfunction as the only dependent variable. 
The vorticity does not enter into this formulation. The algorithm is second order accurate 
in both time and space, it uses a multigrid solver at each time step, and it is extremely 
efficient with respect to use of both CPU time and physical memory. The algorithm is 
extremely robust with respect to Reynolds number, and has been used to directly compute 
incompressible flows with smoothly resolved streamfunction, kinetic energy and vorticity 
contours for Reynolds numbers as high as Re = 100,000 without requiring any subscale 
modelling. Solutions are shown for cavity flows at various Reynolds numbers. 


2. THE ALGORITHM 


In a bounded open region H C R 2 , if 0 is the streamfunction, then the equation for 
two dimensional time dependent viscous incompressible flows with no body force can be 
written as 


3A0 d± A d± <>± A d± 

dt dy dx dx dy 


Re 


A 2 0 = 0, 


( 1 ) 


for x in 0, and t > 0. Notice that this is a single equation for a single scalar unknown, 
and that neither vorticity nor pressure enter into the formulation. The velocity solution is 
obtained as 


u(x,t) 


dtp 

dy' 


and v(x,t) — — 


dtp 

~di' 


( 2 ) 


and it is always divergence free. One of the advantageous features of this formulation is 
that the scalar unknown tp is smoother than the velocity and the vorticity, since they are 
both obtained from the streamfunction by differentiation. The data for the problem in 
this formulation consists of the initial data for tp , and boundary data such as the standard 
set of boundary conditions 


tp(x,t) = /3(x,t), 


and 


dtp 

dr} 


(M) = If (*,«). 


( 3 ) 


for x in dfl. and t > 0, where S- is differentiation in the exterior normal direction at the 

7 7 a rf 

boundary. A steady state form of Equation (1) was previouosly used by Schreiber and 
Keller [7] with path continuation methods for calculating high Reynolds number steady 
cavity flows. 


Let z n be the discrete grid function approximating tp at time t„. Using centered 
spatial differencing throughout, let La approximate the laplacian, let Bi approximate the 
biharmonic operator, and if 


then let Cv approximate the convection terms. We discretize equation (1) as 


Ll(r+1) -^ Bi(r+,) 


_ , At 3 At _ At 

= La(z ) + — Bi(z ) - — Cv(z ) + — Cv(z 1 ). 


( 5 ) 


Notice that equation (5) for z n + 1 is elliptic for all Reynolds numbers, and that it is linear 
with all of the nonlinearities lagged into the source term calculated from data at times 
t n and t n _ i. The diffusion terms are time differenced with a Crank-Nicolson differencing 
scheme, which gives an infinite speed of propogation for the data at each time step through 
the elliptic biharmonic operator, and which does not impose a stability constraint on the 
time step size. The convection terms are time differenced with a lagged second order 
Adams- Bashforth differencing scheme, which does impose the CFL constraint < 1, 

where we assume that Ax = Ay. We have used a CFL number between 0.6 and 0.8 for 
calculations with Reynolds numbers between 100 and 100000, with grid sizes between ^ 
and The local domain of dependance is the large symmetric 13 point discretization 
stencil from the discrete operators Bi and Cv. The velocity components are 

directly recovered using (2) with centered differences, and are both defined at every grid 
point along with the streamfunction approximation .. The discrete velocity solution is 
exactly incompressible, since 6 X (u” .) +S y (w“ .) = 0. This algorithm is related to a primitive 
variable method on a MAC staggered grid by the Finite Difference Galerkin Method (see 
Goodrich and Soh [5]). 


When (5) is used as the finite difference approximation to (l), the problem that must 
be solved at each time step may be written as 

u («" +1 )-^ B i(* n+1 ) = r-"-‘, (6) 

where f n,n_ 1 is the discrete source term from the right hand side of equation (5). We 
have used a banded LU decomposition solver to directly solve equation (6) by back sub- 
stitution at every time step. In order to avoid the large storage overhead required by this 
approach, we currently use a multigrid method at each time step to solve equation (6). 
The multigrid solver factors Bi as two laplacians, it uses a Gauss Seidel or red black Gauss 
Seidel smoothing iteration, a linear restriction and prolongation, and a simple V cycle 
multigrid iteration using three smoothing sweeps while coarsening and none while refining. 
The factoring of Bi as two laplacians follows Linden [6], and introduces u> = Aip only for 
the purpose of having a convergent iteration scheme, but it incidentely produces both 0 
and <jJ as simultaneeous solutions of a coupled Poisson and Laplace’s equation. Note that 
fn.n-i j g a i wa y S calculated from the right hand side of equation (5) using just the discrete 
streamfunction fields z n and z n " 1 . Since w = A0 is only introduced as an intermediate 
variable for solving equation (6) at each time step, its tendency toward a greater sensitivity 
to spatial disturbances and errors is not propagated to the solution 0 at subsequent time 
steps. 
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3. CAVITY CALCULATIONS 


The examples for demonstrating the performance of this algorithm will be some typical 
two dimensional driven cavity calculations. For all driven cavity calculations, irrespective 
of Reynolds number and grid size, between 10 and 15 V cycles are used at each time step 
to reduce the residuals for equation (5) to less than 5 x 10 -12 . 

The performance of the algorithm will be shown by comparision of the calculation of a 
transient square driven cavity flow on an IBM RS/6000 model 530 workstation with three 
different uniform grids. The transient flow is the development from no initial flow with 
an impulsively started lid, from t = 0 to t = 1, with Re = 9600. The three calculations 
are: (a) with a 128 x 128 spatial grid, with 6 grid levels, and with At = requiring 
1.7 MBytes of memory, and 2.5 CPU seconds per time step; (b) with a 192 x 192 spatial 
grid, with 7 grid levels, and with At — requiring 3.8 MBytes of memory, and 5.5 CPU 
seconds per time step; and (c) with a 256 x 256 spatial grid, with 7 grid levels, and with 
At = requiring 6.6 MBytes of memory, and 10.2 CPU seconds per time step. Notice 
first that the CPU time per time step increases linearly with the number of grid points. 
Simulation case (a) requires 1.50 x 10“ 4 seconds per time step per grid point, case (b) 
requires 1.48 x 10~ 4 seconds per time step per grid point, and case (c) requires 1.54 x 10“ 4 
seconds per time step per grid point. The computional requirement for refined grids is 
not increased for each grid point for each time step, but the CFL constraint will require a 
time step size that decreases linearly with the spatial grid size for the refined grid, so that 
more time steps are required for the same nondimensional time interval. If the additional 
time steps that are required for the entire calculation from t = 0 to t = 1 are considered 
as part of the computational cost, then refining the grid from a 128 x 128 to a 192 x 192 
grid actually requires 1.6 times more CPU time per grid point for the entire calculation 
on the finer grid, and refining again from a 192 x 192 to a 256 x 256 grid also requires 
1.6 times more CPU time per grid point for the entire calculation. Notice also that the 
memory requirement increases linearly with the number of grid points. Simulation case 
(a) requires 107.12 Bytes per grid point, case (b) requires 106.97 Bytes per grid point, and 
case (c) requires 104.78 Bytes per grid point. The observed linear increase in CPU and 
memory requirements with an increase in the number of grid points is typical of multigrid 
solvers, since it is well known that the order of computational effort for these solvers is 
0(N ), where N is the total number of grid points. 

The accuracy of the algorithm will be shown by comparing sample results for the 
driven cavity with established results in the literature by Ghia et.al. [1]. Figure (1) presents 
driven cavity results with a 256 x 256 grid for Re — 5000. This data was calculated using 
the time dependent algorithm presented above with At = and is for the state of the 
flow at t = 748.45 when the time evolution has virtually stopped. The cavity has its upper 
lid moving from the left to the right, and all three of its other walls are not moving. The 
length scale is such that the cavity walls have length 1, and the velocity and time scales 
are established by using 1 as the lid velocity. The streamfunction contours in Figure (la) 
show the usual assortment of primary and secondary circulation patterns, but the contours 
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do not show tertiary circulations that are actually resolved in the lower two corners. The 
tertiary recirculating eddy in the lower right hand corner is resolved by 14 x 11 grid points, 
while the smaller tertiary eddy hidden below the streamfunction contours in the lower left 
hand corner is resolved by only 4x3 grid points. Figure (lb) shows the vertical velocity 
component in a horizontal crossection across the middle of the cavity. The calculated 
data is compared with the data in Ghia, et.al. [1], and shows complete agreement with 
the published data. Figure (lc) shows the horizontal velocity component in a vertical 
crossection up the middle of the cavity, and Figure (id) shows the vorticity on the upper 
lid. This data also shows complete agreement with the published data. The acurracy of the 
computed solution is also shown by quantitative comparisons between the calculated and 
published data for the values and locations of the streamfunction minimum and maximum. 
The global streamfunction minimum occurs in the center of the primary circulation, and 
the global streamfunction maximum occurs in the center of the lower right hand secondary 
corner eddy. This data is given in Table 1, along with established data from a calculation 
on a 256 x 256 grid published in Ghia, et.al. [1]. There is excellant agreement between 
our solution and the published data in Ghia, et.al. [1], where both calculations are on a 
256 x 256 grid. Our calculated solution on a 128 x 128 grid also shows very good agreement 
with the published data on a 256 x 256 grid. 

The degree of agreement between our calculated results and the published data for 
the driven cavity at Re = 5000 is also characteristic of our calculations at Re = 7500. For 
Re = 10000 (Goodrich [2]), our calculated results differ dramatically from the published 
steady state solutions for the driven cavity. Instead of converging to a steady state, our 
time dependent calculations show convergence to an unsteady time asymptotic state. In 
fact, our calculations (Goodrich [3]) show a steady state for Re = 8900, and periodic time 
asymptotic solutions with spectra that have a single fundamental frequency for Re — 9000, 
Re = 9500 and Re = 9600. All three of these time asymptotic flows with single fundamental 
frequencies have been calculated on 128 x 128 grids with A t = The periodic flow 
for Re = 9600 has a fundamental frequency of 0.55 ± 0.005 on a 128 x 128 grid. A 
second calculation at Re = 9600 was done on a 192 x 192 grid with A t = ^y, and this 
calculation also produced a periodic time asymptotic flow with one fundamental frequency 
of 0.58 ± 0.005. The data for 3600 < t < 3700 from this calculation on the 192 x 192 grid 
is presented in Figure (2a) as a phase portrait of on the vertical axis versus the total 
kinetic energy on the horizontal axis. This data is for 100 nondimensional time units, or for 
25600 discrete time steps, and for approximately 58 complete periodic cycles of the time 
asymptotic flow. The period of the asymptotic flow state is 1.72 ±0.015, so that each cycle 
requires approximately 440 discrete time steps. Notice the fact that the oval plot appears 
to be a single line. This exact repetition of the phase portrait through 58 cycles shows the 
precise periodicity of the asymptotic state. The qualitative and quantitative agreement of 
the dynamics from these two calculations on separate grids supports the grid independence 
of the periodic time asymptotic solution in the driven cavity at Re = 9600. For higher 
Reynolds numbers we find more complex dynamics. At Re = 9700 on a 128 x 128 grid 
we find two incomensurate fundamental frequencies, and the unsteady time asymptotic 
solution is aperiodic with a discrete spectrum. A phase portrait for this flow is given in 
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Figure (2b), which shows a plot of on the vertical axis versus the total kinetic energy 
on the horizontal axis for 3700 < t < 4500. The data plotted in Figure (2b) is for 800 
nondimensional time units, or for 128600 discrete time steps, and it qualitatively shows 
the dramatically more complex yet highly structered aperiodic time asymptotic state at 
Re = 9700. At Re = 10000 on a 128 x 128 grid, we also find two incomensurate fundamental 
frequencies, and this unsteady time asymptotic solution is also aperiodic with a discrete 
spectrum (see Goodrich [2]). These solutions suggest a first Hopf bifurcation in the square 
driven cavity for 8900 < Re < 9000, and a second one for 9600 < Re < 9700. 

It could very well be argued that all of these unsteady time asymptotic solutions 
in the driven cavity are spurious, since there are established steady solutions for Re = 
10000. On the other hand, the most well known of these solutions (Ghia, et.al. [1], 
Schreiber and Keller [7]) use steady state methods, and unsteady time dependent methods 
for high Reynolds number flows tend to introduce artifical dissipation, which we have 
avoided. Furthermore, Shen [8] has used a Chebychev Tau method to detect a Hopf 
bifurcation in the regularized driven cavity for 10000 < Re < 10500. The use of a spectral 
method with similar results lends plausibility to our claims, since the parabolic velocity 
profile for the regularized cavity lid should transmit less momentum to the cavity interior, 
and the transition to unsteady time asymptotic states could reasonably be expected to 
occur at a higher Reynolds number. There are also periodic solutions in the aspect two 
rectangular driven cavity using our algorithm (Goodrich, et.al. [4]) which show the same 
generic behavior in a deeper cavity and at a lower Reynolds number. Unsteady time 
asymptotic flows are completly consistent with what should be expected in light of the 
current understanding of nonlinear dynamical systems (see [9]). It is possible that our 
algorithm has dispersion or other error characterisitcs that prematurely trigger a Hopf 
bifurcation at too low a Reynolds number, but the cumulative effect of the considerations 
introduced above is that the transition to unsteady time asymptotic states such as we have 
observed should be expected, and that the critical Reynolds numbers for the transitions 
that we have observed are plausible. 

The robustness of the algorithm with respect to Reynolds number will be shown by 
data from a direct computation of a transient square cavity flow with an impulsively started 
lid on a 256 x 256 uniform grid at Re = 25000. Figures (3a-b) show streamfunction and 
vorticity contours at t = 10.25, early in the transient flow development from the initial 
state of rest. The vorticity plot clearly shows a large number of small structures that are 
being created in the boundary layer beneath the wall jet, which is descending from the end 
of the moving lid, and which is separating from the right hand wall. This series of small 
scale structures is ejected from the boundary layer and convected around the large central 
recirculation pattern, to appear either as small recirculating eddies or as waves all around 
the central structure in the cavity. Notice in the vorticity plot that the series of small 
scale structures is wrapped around and folded back into the central structure, resulting in 
a twisted or spiral shaped composite structure that has multiple layers like a fine pastry. 
Figures (3c-d) show streamfunction and vorticity contours at t — 83, and Figures (3e-f) 
show stre amf unction and vorticity contours at t = 84. These plots show the state of the 
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flow after sufficient time has elapsed for the development of all the usual secondary flow 
structures in the cavity. In addition to these secondary structures, we also observe a very 
active collection of tertiary flow structures, that originate in the boundaries and corners 
of the cavity, that are convected and distorted by the local flow while interacting with 
the primary, secondary and tertiary structures, and that ultimately either dissipate or 
combine with other structures. As an example, consider the tertiary eddy that is rising 
along the middle of the left hand wall at t = 84. This structure can be seen at t = 83 as 
the small tertiary structure that is just beginning to pull out of the recirculation complex 
in the lower left hand corner. This tertiary eddy actually originated at about t = 78 in 
the lower right hand corner, and it then separated from that corner to be rolled along the 
lower wall, until it attached to and moved along the upper surface of the lower left hand 
corner complex. For t > 84, this tertiary structure will continue up the left hand wall 
until it attaches to and rolls along the surface of the secondary structure in the upper left 
near the lid, ultimately being absorbed by this secondary eddy. Besides the complexity 
of structures on several scales, at t = 83 and t = 84, we also observe that the vorticity 
in the center of the cavity is not uniform. We see a continuous elastic deformation of the 
primary circulation structure in the center of the cavity, which stretches and moves while it 
rotates. The observed flow clearly does not have a central core of uniformly rotating fluid 
with constant vorticity at t « 83. This computation has not yet been carried further than 
t = 100, so that we do not yet know how the asymptotic flow state can be characterised 
at Re = 25000. The issue of whether or not the 256 x 256 grid is adequate for resolving 
the relevant length scales at Re = 25000 is also not answered by this single transient 
calculation. Nevertheless, we do observe that the streamfunction and vorticity fields are 
smoothly resolved by our algorithm at Re = 25000 with a 256 x 256 grid, and this clearly 
shows the robustness of our algorithm with respect to Reynolds number. 
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The Streamfunction Minimum 


Source 

Grid 

V*m in 

*^min 

I /rn in 

Goodrich 

128 x 128 

-1.15 x 10" 1 

66 _ 132 

128 256 

60 . 
128 

Goodrich 

256 x 256 

-1.18 x 10" 1 

132 

256 

137 

256 

Ghia, et.al. 

256 x 256 

-1.19 x 10- 1 

131 

256 

137 

256 



The Streamfunction 

Maximum 


Source 

Grid 

in 

in 

Vm in 

Goodrich 

128 x 128 

3.44 x 10" 3 

102 _ 204 
128 256 

9 _ 

128 

Goodrich 

256 x 256 

3.13 x 10- 3 

206 

256 

19 

256 

Ghia, et.al. 

256 x 256 

3.08 x 10- 3 

207 

256 

19 

256 


Table 1: STREAMFUNCTION MAX AND MIN 
The square driven cavity at Re=5000 


138 

256 


18 

256 
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Figure 1: Re = 5000, 256 x 256 grid, At = t 
calculated data compared with Ghia, et.al. 


la: Stream Function Contours 


lb: v at y=0.5 



10 





Figure 2: Phase Portraits of V’m 
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2b: Re = 9700, 
128 x 128 grid, 
3700 < t < 4500 
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3c: Streamfunction Contours, 
t=83 


3d: Vorticity Contours, 
t=83 
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